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Abstract 

A thermodynamically consistent damage model for the simulation of progressive 
delamination under variable mode ratio is presented. The model is formulated in the 
context of the Damage Mechanics. The constitutive equation that results from the 
definition of the free energy as a function of a damage variable is used to model the 
initiation and propagation of delamination. A new delamination initiation criterion 
is developed to assure that the formulation can account for changes in the loading 
mode in a thermodynamically consistent way. The formulation proposed accounts 
for crack closure effets avoiding interfacial penetration of two adjacent layers after 
complete decohesion. The model is implemented in a finite element formulation. The 
numerical predictions given by the model are compared with experimental results. 


1 Introduction 


Delamination is one of the most common types of damage in laminated fibre- 
reinforced composites due to their relatively weak interlaminar strengths. De- 
lamination may arise under various circumstances, such as in the case of trans- 
verse concentrated loads caused by low velocity impacts. 

Structural collapse in a composite structure is often caused by the evolution 
of different types of damages created in a local zone of the structure. The 
particular damage modes depend upon loading, lay-up and stacking sequence. 
Delamination is often a significant contributor to the collapse of a structure 



[1], Moore et al. [2] justify the need to account for delamination from the point 
of view of Fracture Mechanics: the energy release rates necessary for failure of 
a composite part from intralaminar and interlaminar damage were computed 
and the results showed that the lowest energy release rate obtained was for 
delamination. 

When other material non-linearities can be neglected, methods based on Lin- 
ear Elastic Fracture Mechanics (LEFM) have been proven to be effective in 
predicting delamination growth. However, LEFM cannot be applied without 
an initial crack. In many situations, stress-based methods have been applied to 
predict the initiation of delaminations, and after the delamination onset, Frac- 
ture Mechanics can be used to describe delamination growth [3] -[4], Techniques 
such as virtual crack closure technique (VCCT) [5] -[9], J-integral method [10], 
virtual crack extension [11] and stiffness derivative [12] have often been used. 
These techniques have in common that the delamination is assumed to prop- 
agate when the associated energy release rate is greater than or equal to a 
critical value [13]. However, difficulties are also encountered when these tech- 
niques are implemented using finite element codes. The calculation of fracture 
parameters, e.g. stress intensity factors or energy release rates, requires nodal 
variable and topological information from the nodes ahead and behind the 
crack front. Such calculations can be done with some effort for a stationary 
crack, but can be extremely difficult when progressive crack propagation is 
involved. 

Another approach to the numerical simulation of the delamination can be 
developed within the framework of Damage Mechanics. Models formulated 
using Damage Mechanics are based on the concept of the cohesive crack model: 
a cohesive damage zone or softening plasticity is developed near the crack 
front. The origin of the cohesive crack model goes back to Dugdale [14] who 
introduced the concept that stresses in the material are limited by the yield 
stress and that a thin plastic is generated in front of the notch. Barenblatt [15] 
introduced cohesive forces on a molecular scale in order to solve the problem 
of equilibrium in elastic bodies with cracks. Hillerborg et al. [16] proposed 
a model similar to the Barenblatt model, but where the concept of tensile 
strength was introduced. Hillerborg’s model allowed for existing cracks to grow 
and, even more importantly, also allowed for the initiation of new cracks. 

Cohesive damage zone models relate tractions to displacement jump at an 
interface where a crack may occur. Damage initiation is related to the inter- 
facial strength, i.e., the maximum traction on the traction-displacement jump 
curve. When the area under the traction/displacement jump curve is equal to 
the fracture toughness, the traction is reduced to zero and new crack surface 
formed. The advantages of cohesive zone models are their simplicity and the 
unification of crack initiation and growth within one model. Moreover, cohe- 
sive zone formulations can also be easily implemented in finite element codes 
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using decohesion elements [IT]- [25] . 


In the formulation of the cohesive models, it is important to control the energy 
dissipation during delamination growth in order to avoid the restoration of 
the cohesive state, i.e. , it is necessary to assure that the model satisfies the 
Clausius-Duhem inequality. There are some models in the literature that can 
be used under mixed-mode conditions [21], [23]- [24], [29]- [34], but some of them 
do not satisfy the Clausius-Duhem inequality under variable-mode loading 
situation. 

A damage model for the simulation of delamination under variable-mode is 
presented in this paper. A new delamination initiation criterion is developed 
from the expression of the critical energy release rate for delamination propa- 
gation under mixed-mode loading presented in [35]. The model is implemented 
into the commercial finite element code ABAQUS by means of a user- written 
decohesion element. 

This paper is structured as follows: first, the formulation of the delamination 
onset and growth model is presented. Afterward, the delamination damage 
model is formulated and the initiation and propagation surfaces are defined. 
The finite element discretization of the boundary value problem is based on 
the variational formulation of the local form of momentum balance. Finally, 
the numerical predictions are compared with experimental results. 


2 Model for delamination onset and propagation 


The boundary value problem, the kinematic equations, and the constitutive 
relations are presented for the formulation of the model for delamination onset 
and delamination propagation. 


2.1 Boundary value problem 


Consider a domain fi, as shown in figure 1(a), containing a crack T c . The part 
of the crack on which a cohesive law is active is denoted by T co /j and is called 
the fracture process zone (FPZ). 

Prescribed tractions, F), are imposed on the boundary Fp, whereas prescribed 
displacements are imposed on r u . The stress held inside the domain, , is 
related to the external loading and the closing tractions rf, tJ in the cohesive 
zone through the equilibrium equations: 

a ijj = 0 on 11 (1) 
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Fig. 1. Body fl crossed by a material discontinuity in the undeformed configu- 
ration 


(7 ijTlj = ti Oil r p 


( 2 ) 


(TijUj = rf = ~T i = (J ijTlj Oil Tcoh 


( 3 ) 


2.2 Kinematics of the interfacial surface 


To develop the necessary kinematic relationships, consider the crack T c shown 
in Figure 1(a) part of a material discontinuity, T^, which divides the domain 
hi into two parts, hl+ and hi- (Figure 1(b)). 

The displacement jump [tq] across the material discontinuity can be writ- 
ten as: 

Ikl = u i - u ~i ( 4 ) 

where uf denotes the displacement of the points on the surface of the material 
discontinuity of the part hl± of the domain. 

The fundamental problem introduced by the interfacial surface is how to 
express the virtual displacement jumps associated to the surfaces T d ± as a 
function of the virtual displacements. Consider a three-dimensional space with 
Cartesian coordinates X*, i = 1,2,3 . Let the Cartesian coordinates xf de- 
scribe the deformation of the upper and lower surfaces r^± in the deformed 
configuration. Any material point on r^± in the deformed configuration is 
related to its undeformed configuration through: 

xf = Xi + uf (5) 

where uf are displacements quantities with respect to the fixed Cartesian 
coordinate system. Then, the coordinates Xi of the midsurface can be written 
as: 

Xi = Xi + - (uf + uf ) (6) 
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Undeformed configuration Deformed configuration 


Fig. 2. Interfacial surface deformation 

The components of the displacement jump vector are evaluated at the mid- 
surface T d , which is coincident with T d in the undeformed configuration (see 
figure 2). The midsurface coordinate gradients define the components of the 
two vectors, v v . and v^., that define the tangential plane at a given point, P: 


V V i — X i,V 

(7) 

V ti = 

(8) 


where r], £ are curvilinear coordinates on the surface T d . Although v v and 
are generally not orthogonal to each other, their vector product defines a 
surface normal. Therefore, the local normal coordinate vector is obtained as: 

V„ = V£ x V,, ||v € x vj _1 (9) 

The tangential coordinates are then obtained as: 

v - = v dl v dl _1 ( 10 ) 

v t = v n x v s (11) 

The components of v n , v s and v t represent the direction cosines of the local 
coordinate system in the global coordinate system at a material point P G T d . 
The director cosines define the following rotation tensor 0 m *, expressed in 
Voigt notation as: 

© = [v n ,v s ,Vt] (12) 

0 is orthogonal and relates the local coordinate system to the fixed coordinate 
system. Using the rotation tensor, the normal and tangential components of 
the displacement jump tensor expressed in terms of the displacement held in 
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global coordinates are: 


(13) 


A m — ©mi J^i] 

A m is the displacement jump tensor in the local coordinate system. 

The deformed differential surface area of the midsurface cit'd is expressed in 
the form: 

clt d = ||v c x vj clT° d = JdT° d (14) 

where clF d is the differential undeformed surface area. 


2.3 Constitutive laws 

A constitutive law relating the cohesive tractions, Tj , to the displacement 
jump in the local coordinates, A, , is required for modelling the constitutive 
behaviour of the material discontinuity. The constitutive laws in the material 
discontinuity may be formally written: 

Tj = T (A;) (15) 

tj = Dp A, (16) 

where Djf n is the constitutive tangent stiffness tensor. 

A new constitutive model relating the displacement jumps to the tractions, 
and based on Damage Mechanics is proposed. 

The delamination model proposed follows the general formulation of Contin- 
uum Damage Models proposed by Simo and Ju [36]- [37] and Mazars [38]. 

The free energy per unit volume of the interface is defined as: 

V>(A,d) = (l-d)^°(A) (17) 

where d is the damage variable and is a convex function in the displacement 
jump space defined as: 

V-° (A) = t A,Dl A, (18) 

Negative values of A 3 do not have any physical meaning because interpene- 
tration is prevented by contact. Therefore negative values of A 3 should not 
have any influence in the variation of the free energy of the interface. Thus, a 
modification of equation (17) is proposed to prevent interfacial penetration of 
the two adjacent layers after complete decohesion. The expression for the free 
energy per unit volume proposed is: 

ip (A, d) = (1 - d) V>° (A*) - d-0 0 ( 5 3i (-A 3 )) (19) 
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where (•) is the MacAuley bracket defined as (x) = | (x + |x|) and Sij is the 
Kronecker delta. The constitutive equation for the interface is obtained by 
differentiating the free energy with the displacement jumps: 


A = = (1 - d) D[ ;.A, - d D^ 3j (-As) (20) 

In the present model, the initial stiffness tensor, , is dehned as: 

D» = S,,K (21) 

where the scalar parameter K is a penalty stiffness. In the expanded form, the 
constitutive equation can be written in Voigt notation as: 


\ 

Tl 


r \ 

Ar 


0 

r 2 

> = (1 - d) K < 

A2 

> — dA" < 

0 

^3 


A3 

\ / 


<-A 3 ) 


(22) 


The energy dissipation during damage evolution, H , represented in Figure 3 
for single-mode loading, can be obtained from: 


d'ip ■ 

= -d > 0 

dd ~ 


(23) 


The model dehned by equation (20) is fully determined if the value of the 



Fig. 3. Energy dissipation during damage evolution 

damage variable d can be evaluated at every time step of the deformation 
process. For that purpose, it is necessary to define a suitable norm of the 
displacement jump tensor, a damage criterion, and a damage evolution law, 
as will be described in the following sections. 


7 



2.3.1 Norm of the displacement jump tensor 

The norm of the displacement jump tensor is denoted as A and is also called 
equivalent displacement jump norm. It is used to compare different stages of 
the displacement jump state so that it is possible to define such concepts as 
‘loading’, ‘unloading’ and ‘reloading’. The equivalent displacement jump is a 
non-negative and continuous function, defined as: 

A = \J (A 3) 2 + (A shear) 2 (24) 

where A 3 is the displacement jump in mode I, i.e. , normal to midplane, and 
A s hear is the euclidean norm of the displacement jump in mode II and in mode 
III: 

A shear = \] (Ai)" + (A 2 )" (25) 


2.3.2 Damage criterion 

The damage criterion is formulated in the displacement jump space. The form 
of this criterion is: 


F (A*, r*) := A* - r* < 0 Mt > 0 (26) 

where t indicates the actual time and r f is the damage threshold for the current 
time. If r° denotes the initial damage threshold, then F > r° at every point 
in time. Damage initiation is produced when the displacement jump norm, A, 
exceeds the initial damage threshold, r°, which is a material property. 

A fully equivalent expression for equation (26) that is more convenient for 
algorithmic treatment is [39]: 

F (A*, r*) := G (A*) - G (V) <0 Mt > 0 (27) 

where G(-) is a suitable monotonic scalar function ranging from 0 to 1 . G(-) 
will define the evolution of the damage value, and will be presented in the 
next section. 


2.3.3 Damage evolution law 

The evolution laws for the damage threshold and the damage variable must be 
defined in the damage model. These laws are defined by the rate expressions: 

r = /i (28) 


d 


. dF (A, r) 

OX 


. dG (A) 
fl ^x~ 


(29) 
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where /i is a damage consistency parameter used to define loading/unloading 
conditions according to the Kuhn- Tucker relations: 

/i> 0 ; F(A*y)<0 ; /IF (A*, r*) = 0 (30) 

From the previous equations it is easy to prove [36] that the evolution of the 
internal variables may be explicitly integrated to render: 

r* = max |r°, max A s | 0 < s < t (31) 


d* = 



(32) 


which fully describes evolution of the internal variables for any loading/ un- 
loading/ reloading situation. The scalar function G (•) defines the evolution 
of the damage value. For a given mixed-mode ratio, (3, the function proposed 
here is defined as, 


G (A) = 


A f (A - A 0 ) 


(33) 


A (A / - A 0 ) 

Equation (33) defines the damage evolution law by means of a bilinear consti- 
tutive equation (see Figure 4), where A 0 is the onset displacement jump and 
it is equal to the initial damage threshold r°. The initial damage threshold 
is obtained from the formulation of the initial damage surface or initial dam- 
age criterion. A? is the final displacement jump and it is obtained from the 
formulation of the propagation surface or propagation criterion. 


Onset 



Fig. 4. A bilinear constitutive equation for the decohesion element for a mixed mode 
loading situation 


It is therefore necessary to establish the delamination onset and propagation 
surfaces for the complete definition of the damage model. Delamination on- 
set and propagation surfaces and the damage evolution law fully define the 
constitutive equations. 
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The constitutive equations for the interfacial surface are normally developed in 
a phenomenological way, i.e. , satisfying empirical relations that are obtained 
using experimental results. There are several types of constitutive equations 
used in decohesion elements: Tvergaard and Hutchinson [40] purpose a trape- 
zoidal law, Cui an Wisnom [41] a perfectly plastic rule, Needleman first pre- 
sented a polynomial law, [27], and later an exponential law [28]. Goyal et al. [42] 
adopted Needleman’s exponential law to account for load reversal without ma- 
terial healing. 


In this paper the law used within the decohesion elements is a bilinear law 
[21], [24], [43]. The bilinear law is the most commonly used cohesive law due to 
its simplicity. One drawback of the bilinear law is that the traction-displacement 
jump relation is discontinuous at its maximum. The discontinuities in the 
traction-displacement jump relation can be avoided using continuous func- 
tions. However, even for such continuous functions, the discontinuity is un- 
avoidable when modeling loading-unloading cycles. 


For a given mixed-mode ratio, /?, defined as: 


P 


A 


shear 


^ shear H - (^3) 


(34) 


the bilinear constitutive equation is defined by a penalty parameter, K, the 
damage value, d, the mixed- mode damage initiation, A 0 and the total decohe- 
sion parameter, . These last two values are given by the formulation of the 
onset and the propagation criterion which takes into account the interaction 
between different modes and their value depends on the mixed-mode ratio (3. 
The penalty parameter K assures a stiff connection between two neighboring 
layers before delamination initiation. The penalty parameter should be large 
enough to provide a reasonable stiffness but small enough to avoid numerical 
problems, such as spurious tractions oscillations [44], in a finite element anal- 
ysis. Daudeville et al. [45] have proposed the definition of the penalty stiffness 
as a function of the interface thickness, t, and elastic moduli of the interface. 
Zou et al. [4], [46] have proposed a value for the penalty parameter between 10 5 
and 10 7 times the value of the interfacial strength per unit length. Camanho 
and Davila [24], [47] successfully used a value of 10 6 N/mm 3 , which is smaller 
than the values proposed by Daudeville et al. [45] and Zou et al. [4], [46]. 


Propagation criterion 

The criterion used to predict delamination propagation under mixed-mode 
loading conditions is usually established in terms of the energy release rate 
and fracture toughness. It is assumed that when the energy release rate, G, 
exceeds a critical value, the critical energy release rate G c , delamination will 
grow. 
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The most widely used criterion to predict delamination propagation under 
mixed-mode loading, the ” power law criterion” is normally established in 
terms of a linear or quadratic interaction between the energy release rates 
[48]. However, Cama nh o et ah [24] shown that using the expression proposed 
by Benzeggagh and Kenane [35] for the critical energy release rate for a mixed- 
mode ratio is more accurate for epoxy composites. Using the expression by 
Benzeggagh and Kenane in [35], the propagation criterion can be written as: 

G = G Ic + (G IIc - Gf c ) (35) 

where G = Gi + G s h ear is the energy release rate under mixed-mode loading 
and G s hear = G n + Gm is the energy release rate for shear loading proposed 
by Li [49], [50]. 

The propagation surface in the displacement jump space is defined through 
the final displacements, which are defined from the pure mode fracture tough- 
ness ( Gic , Guc , Gnic ) and considering that the area under the traction- 
displacement jump curves is equal to the corresponding fracture toughness, 
i.e.: 

G c = ^KA°A f (36) 

For a given mixed-mode ratio, /?, the energy release rates corresponding to 
total decohesion are obtained from: 

G, = Ua» 03) A( (ft) (37) 


G, hmr = l -K 03 ) (/?) (38) 

where A° shear (f3) and (/?) are respectively the shear and normal displace- 
ment jump corresponding to the onset of softening under mixed- mode loading, 
and A{ hear (/ 3 ) and A 3 (/ 3 ) are respectively the shear and normal displacement 
jump corresponding to the total decohesion under mixed-mode loading. 


From (34): 

AL„ r 08) = A" 03 ) -^L- (39) 

AL„ ((3) = A' (0) LL- (40) 

Using equations (39) and (40) in (37) and (38) the ratio between G ^y ar can 
be established in terms of (3: 


r shear 


P 2 

1 + 2(3 2 - 2/3 


(41) 
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Since the ratio G ^y ar is only a function of the mixed-mode ratio /?, hencefor- 
ward this ratio is named as B: 


B 


r shear 

Gt 


(42) 


Using equation (36), (41) and (42) in equation (35) the propagation criterion 
is obtained in the displacement jump space as: 


A f 


+ (A1.AL- - A§A0 [BY 1 

A° 


(43) 


Initial damage surface 

Under pure mode I, mode II or mode III loading, delamination onset occurs 
when the corresponding interlaminar traction exceeds its respective maximum 
interfacial strength, Under mixed-mode loading, an interaction be- 

tween modes must be taken into account. Few models take into account the 
interaction of the traction components in the prediction of damage onset. The 
models that account for the interaction of the traction components are usually 
based on Ye’s criterion [51], using a quadratic interaction between modes: 



However, experimental data for the initiation of delamination under mixed- 
mode is not readily available and consequently, failure criterion that can pre- 
dict the initiation have not been fully validated. 

The criterion for propagation is often formulated independently of the crite- 
rion for initiation. In this paper, a link between propagation and initiation 
is proposed. Since delamination is a fracture process, the initiation criterion 
proposed in this paper evolves from the propagation criterion and the damage 
evolution law. The isodamage surface for a damage value equal to 1 corre- 
sponds to the propagation surface obtained from equation (35). Then, the 
isodamage surface for a damage value equal to 0 is the initial damage surface. 
With these assumptions, the criterion for delamination initiation proposed 
here is: 

( T °) 2 = ( r 3) 2 + (ri) 2 + (r 2 ) 2 = (rlf + ([T° hear f - (rg) 2 ) [By (45) 
In the displacement jump space, the criterion becomes: 

(A 0 ) 2 = (A 3) 2 + (A ,) 2 + (A 2 ) 2 = (A ») 2 + ((Al „) 2 - (A 3 ) 2 ) [BY (46) 
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The initiation criterion developed here and summarized by equation (45) is 
compared with Ye’s criterion and with a maximum traction criterion, which 
do not take into account mode interaction. The surfaces obtained by the dif- 
ferent criterions are represented in Figure 5. The values predicted by the new 
criterion are very close to Ye’s criterion, that has been successfully used in 
previous investigations [24], 



Maximum traction 
Ye’s criterion 
Proposed criterion 


Fig. 5. Comparision between Ye’s criterion, a maximum traction criterion and the 
new proposed criterion 


The formulation presented assures a smooth transition for all mixed-mode ra- 
tios between the initial damage surface to the propagation surface through 
damage evolution. In Figure 6 it is represented the evolution of the damage 
surface from the damage initiation surface to the propagation surface for pos- 
itive values of displacement jumps. 


2-4 Formulation of the constitutive tangent tensor 


According to equation (16), the constitutive tangent tensor needs to be defined 
for the formulation of the constitutive laws in the material discontinuity and 
for the later numerical implementation of the proposed model. The constitu- 
tive tangent tensor is obtained from the differentiation of the secant equation 
( 20 ): 


t i DijAj 5 i;j K 


1 + 5. 


3 j- 


b A >) 


j j 


Ajd 


(47) 
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Fig. 6. Damage evolution surface in the relative displacement's space 


where is defined as: 


Dij f}j d\ 


1 - d 1 + 5; 


3 j- 


a* . 


(48) 


The evolution of the damage variable d only occurs for loading situations. 
Then, the evolution of the damage variable can be written as: 


d = 


G(A) = 5ffi>A ,r<A<A ' 

0 , r > A or Al < A 


(49) 


where the variation of the function G is obtained assuming that the variation 
of the final displacement jump A-f and the onset displacement jump A 0 with 
the mixed-mode ratio B are not significant for the time increment taken: 


<9G (A) _ A-f A 0 1 
dX ~ A/ - A° A 5 


(50) 


The evolution of the displacement norm is obtained from equation (24): 


A 



A k 

T 



A fc J 




(51) 


Using equations (49) through (51), equation (47) can be written as: 


Ti 


D tan A • 


(52) 


14 


7~)tan i 

u ij ~ ' 


{ Dij - K [l + 5. 


(=M 

3 i A, 


D, 


where H is a scalar value given by: 


l + S 3i tM.]HA l A J } ,r < X < A f 

, r > X or < A 
(53) 


H = 


At A 0 1 
Af - A° F 


(54) 


3 Finite element discretization - computational model 


To transform the strong form of the boundary value problem into a weak form 
better suited for finite element computations, the displacement U{ must belong 
to the set U of the kinematically admissible displacement fields which allows 
for discontinuous displacements across the boundary T^ of the delamination. 
The weak form of the equilibrium equation is [52]: 

^ ( TijSeijdQ + j Ti5A m JclT ° = J FiSvidT + j pbAvidVl \/vi G U (55) 

where bj, are the body forces, Fi are surface forces, and equation (14) has 
been used to relate the initial and the deformed midsurface of the material 
discontinuity. 


The discretization of the domain has been performed by the discretization of 
the whole domain with standard volume elements. However, the surfaces 
surrounding the potential delamination T,/ are discretized with decohesion el- 
ements [24], The discretized formulation is divided in the two domains consid- 
ering no formal coupling between the continuous and the discontinuous parts 
of the deformation in the expression for the free energy of the interface [53]. 
Since the decohesion elements used have zero thickness, the body forces are 
neglected in these elements. Moreover, it is assumed that no external surface 
loads are applied. Therefore, the weak form of the internal virtual work for 
the elements in the surfaces surrounding the potential delaminations can be 
written as: 


r d 


r^A m J d r|J - 0 


(56) 


The displacements and displacement gradients for the decohesion elements are 
approximated as: 

«i|n, = KC (57) 

M In. = N' K q‘ K , (58) 
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where 



Nfc *'er i 


(59) 


*'er, 

where, q^i is the displacement in the i direction of the K node of the el- 
ement, Nk are standard Lagrangian shape functions [52], N e R are standard 
Lagrangian shape functions dehned for the decohesion elements [24], 


Using equation (58) in equation (56) the weak form of the internal virtual 
work at current time t E R + , is given by: 

[ r iB imK 8q K i J dT° d = 0 (60) 

where B m iK is the discrete operator that relates the displacement jump in 
local coordinates to nodal displacements, 

BimK&QKi (61) 

The system of equations given by equation (60) forms the basis for the assumed 
displacement finite element procedure. 


3.1 Discretization of the interfacial surface 


3.1.1 Element kinematics 


According to equation (57), the displacement held, tq, and the undeformed 
material coordinate, A*, associated to the surfaces T d + are interpolated as 
follows: 

uf = Nk'Ik. (62) 


Xf = "KpKi (63) 

where q^ t are the nodal displacement vector and p Ri are the undeformed 
material nodal coordinate vector. Note that the values of p Ri and p Ri can be 
different in the case that an initial crack exists. Using these equations, the 
material coordinates of the interfacial midsurface are: 


Xi = -N Ki (pi, + p Ki + qi, + q K (j 


(64) 


The components of the two vectors that define the tangential plane can be 
written as: ^ 

% = x i>ri = N KijV - (pfti + p Kl + + fe) (65) 


( PKi + PKi + ( lh + Qki) 


( 66 ) 
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Using (58) and (12), the displacement jump can then be obtained in local 
coordinates as: 

A m = ®im,NkQKi = (67) 


3.1.2 Element internal force vector and tangent stiffness matrix 


The internal force vector of the interface element obtained from equation (60) 
is given by: 


ffi 



T i-BimK J dT c 


( 68 ) 


where -B* m A can be written as: 


li irn l\ 


SAr, 


dB 


jmP 


dqxi 9q K i 


q P j + B imK 


(69) 


The softening nature of the decohesion element constitutive equation causes 
difficulties in obtaining a converged solution for the non-linear problem. Nu- 
merical algorithms such as Newton-Raphson method are used to solve the 
nonlinear problem. Therefore, the tangent stiffness matrix must be defined. 
The tangent stiffness matrix stems from the linearization of the internal force 
vector and it is obtained using Taylor’s series expansion about the approxima- 
tion q^i [25]. The tangent stiffness matrix, K KlZv , for the decohesion element 
is: 


K 


KiZv 


dfia 


dr j 


dq Zv -Jr d dq Zv 


BimK J dT C 


dB, 


imK 


dqzv 


TiJdT* + 


r dJ 
r d dq Zv 


T iBi m KdT' d 


( 7 °) 

The computation of the tangent stiffness matrix is intensive and a very ac- 
curate expression is not required. Therefore, the second and the third term 
in the R.H.S of equation (70) are neglected. Thus, the approximate tangent 
stiffness matrix is: 


K KlZv « / p-B imK J<n* d = [ B vjZ D^B imK JdT° d (71) 
■n ,i dq Zv J r d J 

where is the material tangent stiffness matrix, or constitutive tangent 
tensor used to define the tangent stiffness matrix. The constitutive tangent 
tensor is defined in 2.4. 


4 Comparison with experimental studies 


The formulation proposed here was implemented in the ABAQUS Finite El- 
ement code [54] as a user- written element subroutine (UEL). To verify the 
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element under different loading conditions, the double cantilever beam (DCB) 
test, the end notched flexure (ENF) test, and mixed- mode bending (MMB) 
tests are simulated. The numerical predictions are compared with experimen- 
tal data. The DCB test consists of pure mode I delamination. The ENF tests 
measure pure mode II interlaminar fracture toughness, and the MMB de- 
laminate under Mixed mode I and II. In the absence of mode III loading, 
G s hear = Gjj. To investigate the accuracy of the formulation in the simu- 
lation of delamination DCB, ENF and MMB simulations are conducted for 
PEEK/APC2, a thermoplastic matrix composite material. 


4-1 Mode I, mode II and mixed-mode I and II delamination growth for a 
PEEK composite 


The most widely used specimen for mixed-mode fracture is the mixed-mode 
bending (MMB) specimen shown in Figure 7, which was proposed by Reeder 
and Crews [55] , [56] and later re-designed to minimize geometric nonlinearities 
[57]. This test method was recently standardized by the American Society for 
Testing and Materials [58]. 



Fig. 7. MMB test specimen 


The main advantages of the MMB test method are the possibility of using 
virtually the same specimen configuration as for mode I tests, and the ca- 
pability of obtaining different mixed-mode ratios, ranging from pure mode 
I to pure mode II, by changing the length c of the loading lever shown in 
Figure 7. An 8-node decohesion element is used to simulate DCB, ENF and 
MMB tests in unidirectional AS4/PEEK carbon-fiber reinforced composite. 
The specimens simulated are 102-mm-long, 25.4-mm-wide, with two 1.56-mm- 
thick arms. The material properties are shown in Table 1, and a penalty stiff- 
ness K = 10 6 N/mm 3 is used. 
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Table 1. Properties for PEEK/ASf 


En 

E22 — E33 

G12 = G13 

g 23 

Vl 2 = vvz 

122.7 GPa 

10.1 GPa 

5.5 GPa 

3.7 GPa 

0.25 

^23 

Gic 

Guc 

s 

coo 

r° - r° 

0.45 

0.969 kJ/m 2 

1.719 kJ/m 2 

80 MPa 

100 MPa 


The experimental tests were performed at different ratios, ranging from 
pure mode I loading to pure mode II loading. The initial delamination length 
of the specimens (oo) and the mixed- mode fracture toughness obtained exper- 
imentally are shown in Table 2. 


Table 2. Experimental data 


Gu/Gt 

0% (DCB) 

20% 

50% 

80% 

100% (ENF) 

G c [ kJ/m 2 ] 

0.969 

1.103 

1.131 

1.376 

1.719 

o 0 [mm] 

32.9 

33.7 

34.1 

31.4 

39.2 


Models using 150 decohesion elements along the length of the specimens, and 
4 decohesion elements along the width, are created to simulate the ENF and 
MMB test cases. The initial size of the delamination is simulated by plac- 
ing open decohesion elements along the length corresponding to the initial 
delamination of each specimen (see Table 2). These elements are capable of 
dealing with the contact conditions occurring for mode II or mixed-mode I 
and II loading, therefore avoiding interpenetration of the delamination faces. 
The model of the DCB test specimen uses 102 decohesion elements along the 
length of the specimen. The different GII/GT ratios are simulated by apply- 
ing different loads at the middle and at the end of the test specimen. The 
determination of the middle and end loads for each mode ratio is presented in 
[24] . The experimental results relate the load to the displacement of the point 
of application of the load P in the lever (load-point displacement, Figure 7). 
Since the lever is not simulated, it is necessary to determine the load-point 
displacement from the displacement at the end and at the middle of the spec- 
imen, using the procedure described in [24], The B-K parameter r] = 2.284 
is calculated by applying the least-squares fit procedure proposed in [24] to 
the experimental data shown in Table 2. Figure 8 shows the numerical predic- 
tions and the experimental data for all the test cases simulated, and Table 3 
shows the comparison between the predicted and experimentally determined 
maximum loads. It can be concluded that a good agreement between the 
numerical predictions and the experimental results is obtained. The largest 
difference (—8.1%) corresponds to the case of an MMB test specimen with 
jr 1 - = 20%. This fact is not surprising, since the largest difference between the 
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800 





Fig. 8. Numerical and experimental results 

fracture toughness experimentally measured and the one predicted using the 
B-K criterion occurs for < jf J - = 20% (see [24]). 


Table 3 Comparison of the maximum load 


Gh/Gt 

Predicted [N] 

Experimental [N] 

Error (%) 

0% (DCB) 

152.4 

147.5 

3.4 

20% 

99.3 

108.1 

-8.1 

50% 

263.9 

275.3 

-4.2 

80% 

496.9 

518.7 

-4.2 

100% (ENF) 

697.1 

748.4 

-6.9 


5 Concluding remarks 


A thermodynamically consistent damage model for the simulation of progres- 
sive delamination based on Damage Mechanics was presented. A constitutive 
equation for the interface was derived from the variation of the free energy of 
the interface. The resulting damage model simulates delamination onset and 


20 


delamination propagation. The constitutive equation proposed uses a single 
variable to track the damage at the interface under general loading condi- 
tions. A new initiation criterion that evolves from the Benzeggagh-Kenane 
propagation criterion has been developed to assure that the model accounts 
for changes in the loading mode in a thermodynamically consistent way and 
avoids material healing. The damage model was implemented in the finite el- 
ement code ABAQUS by means of a user- written decohesion element. The 
material properties required to define the element constitutive equations are 
the interlaminar fracture toughnesses, the penalty stiffness, and the strengths 
of the interface. In addition, a material parameter rj, which is determined from 
standard delamination tests, is required for the Benzeggagh-Kenane mode in- 
teraction law. 

Three examples were presented that test the accuracy of the method. Simula- 
tions of the DCB and ENF tests represent cases of single-mode delamination. 
MMB tests were simulated at various proportions of mode I and II loading 
conditions. The examples analyzed are in good agreement with the test results 
and they indicate that the proposed formulation can predict the strength of 
composite structures that exhibit progressive delamination. Although the ex- 
amples presented in this work were obtained for composite specimens contain- 
ing pre-existing delaminations, the formulations can be extended to composite 
structures without any pre-existing defects. 


APPENDIX A 
A Algorithm 


In this section the algorithm implementation of the previous model is outlined. 
The algorithm contemplates two different options for the computation of the 
initiation and propagation surfaces. The first option uses the expression for 
the critical energy release rate formulated by Benzeggagh and Kenane [35]. 
The second option uses a Power-law form for the failure criterion, which is 
widely used in the literature for the simulation of delamination. The Power- 
law criterion is given by: 


Gj 


G 


Ic 


GuJ 


(A.l) 


Initial data for time t+1 

Material properties: G IC , G IIC , G UI c , E, r/, r° shear , 
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Current values: A 3 , A shear, d* 

(1) Determine mixed-mode ratios 


P = 


*shear 


B = 


(A 3 ) -|- A shear 

P 2 

1 + 2 (3 2 - 2(3 

(2) Determine pure mode onset displacements 

T° 

A° = -^ i = 3, shear 
1 K 

(3) Determine mixed-mode onset displacement 

Benzeggagh — Kenane (BK) 


A° = ^(A §)' + [(A ° hear y - (A§) J [B] v 

Power law (PL) 

A 0 = \/l+2/3 2 — 2/3A?j A° fapnr 


(d-/3) A ^ 


\ 277 


■(^s) 


2 V 


(4) Determine mixed-mode final displacements 


BK 

A' = [G Ic + (Gnc ~ G Ic ) [B] 11 } 

PL 


A / = 


2 [l-|-2/3 2 — 2/3 


K A° 


( W + /£L)’ 

V Wc / \ Gr 1 1 c / 


(5) Evaluate displacement jump norm 

A=\/ (ZA 3 )" + (A s ^ ear )‘ 

( 6 ) Update internal variables 

A°A f 


r f = 


r t+1 = max 


A f - d* [A / - A 0 ] 
{r‘,A«} 

,1+1 A / (r 1+1 - A°) 

r i+i (A / - A") 


(A.2) 

(A.3) 

(A.4) 


(A.5) 


(A. 6 ) 


(A.7) 

(A. 8 ) 

(A.9) 
(A. 10) 
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(7) Compute tangent stiffness tensor 


D 


tan 


D 


tan 


+ D 


tan 

k 


where 


D 


tan 


(1 — d) k 0 0 

0 (1 - d) k 0 

0 0 (l-d^)fc 


and Dj, an depends on the loading/unloading conditions: 
IF (A' +1 < r* or A I < A ,+1 ) 


(A-11) 


(A- 12) 




Bf n = 0 


(A.13) 

IF ( 

V < A m < Af) 






-kHA^ 

-kHA 1 A 2 

-kHA.As | 

(As)\~ 
t a 3 ) 


B f k an = 

-kHA 2 A 1 

—kHA 2 A 2 

—kHA 2 Az | 

Qa 3 >\ 
t a 3 ) 



-kH A 3 A! -kH A s A 2 -kH A 3 A 3 (1M) 

(A. 14) 


where for a bilinear constitutive law H is given by equation 54. 
(8) Compute tractions 




' \ 
Ar 


0 

T2 

1 

A2 

> — dA; < 

0 

r 3 


A3 

k > 


(-As) 


(A-15) 
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